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ABSTRACT 

Galactic cosmic ray (CR) acceleration to the knee in the spectrum at a few PeV is only 
possible if the magnetic field ahead of a supernova remnant (SNR) shock is strongly 
amplified by CR escaping the SNR. A model formulated in terms of the electric charge 
carried by escaping CR predicts the maximum CR energy and the energy spectrum of 
CR released into the surrounding medium. We find that historical SNR such as Cas 
A, Tycho and Kepler may be expanding too slowly to accelerate CR to the knee at 
the present time. 

Key words: cosmic rays, acceleration of particles, shock waves, magnetic field, ISM: 
supernova remnants 
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1 INTRODUCTION 

During diffusive shock acceleration (DSA) cosmic rays (CR) 
gain energy by repeatedly passing back and forth between 
the upstream and downstream plasmas (Krymskii 1977, Ax- 
ford et al 1977, Bell 1978, Blandford & Ostriker 1978). CR 
diffuse ahead of the shock to form a precursor with an expo- 
nential scaleheight D u /u s where u„ is the shock velocity and 
D u is the CR diffusion coefficient upstream of the shock. The 
average dwell-time spent upstream of the shock between suc- 
cessive shock crossings is 4D u /cu s for relativistic particles 
(Bell 2012). This, along with the corresponding downstream 
dwell-time, determines the rate at which CR are accelerated. 
Lagage & Cesarsky (1983a,b) showed that the time taken for 
CR acceleration is t acce i = 4D u /Ug+4Dd/Ud where Dd is the 
downstream diffusion coefficient and Ud is the downstream 
fluid velocity in the shock rest frame. Since Ud = u s /4 for 
a strong shock it might appear that the downstream dwell- 
time determines the acceleration rate, but we can expect 
Dd -C D u partly because the magnetic field is increased by 
compression by the shock, partly because the compressed 
downstream magnetic field is more closely perpendicular to 
the shock normal, and partly because the downstream field 
is disturbed and more irregular after passing through the 
shock. If the downstream and upstream dwell-times are the 
same t aC cei = &D u /u 2 s . A further common assumption is that 
Bohm diffusion applies: D u — cr g or D u — cr 9 /3 where r g 
is the CR Larmor radius. This assumes a diffusion model 
in which CR are scattered by irregularities in the magnetic 
field such that the scattering mean free path is of the order of 
the CR Larmor radius. There is some observational evidence 
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for Bohm diffusion (Stage et al 2006, Uchiyama et al 2007). 
Furthermore if the mean free path were much larger than a 
Larmor radius acceleration by SNR to the knee in the Galac- 
tic CR spectrum would be very difficult. The maximum CR 
energy is determined by the condition that t aC cei cannot ex- 
ceed the age tsNR of an SNR (Lagage & Cesarsky 1983a,b). 
Assuming D u = cr g , the maximum CR energy T ma x in eV 
is given by T max = 0.03B M GW?iioooPeV where B^q is the 
upstream magnetic field in microGauss, 117 is the shock ve- 
locity in units of 10,000 km s~ x and £1000 is the SNR age 
in 1000's of years. For a typical young SNR expanding into 
the interstellar medium without magnetic field amplification 
B M G = 3, m = 0.5 and t 1000 = 0.4, giving T max = O.OlPeV 
which falls a factor of ~ 100 short of that required to explain 
the Galactic CR spectrum. SNR in the Sedov phase do not 
fare better. Their shock velocities decrease in proportion to 
time -3 / 5 and radius -3 / 2 , so little benefit accrues from their 
larger age and radius. This posed a serious problem for dif- 
fusive shock acceleration as an explanation of the Galactic 
spectrum until it was shown that a plasma instability driven 
by streaming CR in the upstream precursor could amplify 
the magnetic field ahead of the shock and facilitate rapid ac- 
celeration to higher energies (Lucek & Bell 2000, Bell 2004, 
2005). 

The phenomenon of magnetic field amplification pro- 
vides a mechanism by which the CR energy can be raised 
significantly beyond O.OlPeV, but there remains the ques- 
tion of why the fields are amplified to the observed magni- 
tude, up to 100's /j,G in the historical SNR (Vink & Laming 
2003, Berezhko et al 2003, Volk et al 2005), and why CR are 
accelerated to a few PeV rather than 0.1 or 10 PeV. We try 
to answer these questions by examining the self-consistent 
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interaction between streaming CR, behaving kinetically, and 
the upstream plasma behaving magnetohydrodynamically. 

It has been known for many years that the escape of CR 
upstream of the shock is an important part of the overall ac- 
celeration process as discussed below in the final paragraphs 
of section 3. We find that the combined CR-MHD system 
organises itself to allow a suitable number of CR to escape 
upstream. The CR drive magnetic field amplification which 
in turn regulates the number of escaping CR. If a smaller 
number of CR escaped, the magnetic field would be insuffi- 
ciently amplified to confine and accelerate the CR. If a larger 
number of CR escaped, the magnetic field would grow too 
rapidly to allow their escape. Hence a self-regulating system 
is set up that determines the number and maximum energy 
of escaping CR. 

The paper is organised as follows. Sections 2 and 3 
present approximate calculations showing how a limit on 
the CR energy is placed by the need for CR to drive mag- 
netic field amplification by escaping upstream. Sections 4 
to 7 describe Vlasov-Fokker-Planck (VFP) simulations that 
support the arguments of sections 2 and 3. Sections 8 to II 
apply the results to supernova remnants and the Galactic 
cosmic ray spectrum. Readers unfamiliar with VFP simula- 
tions may wish to read sections 1 to 3 and 8 to 11 before 
returning to the computational validation and illustration in 
sections 4 to 7. 



2 CONDITIONS FOR STRONG MAGNETIC 
FIELD AMPLIFICATION 

We assume that magnetic field is generated by the non- 
resonant hybrid (NRH) instability described by Bell (2004). 
This is one of a class of plasma instabilities driven by CR 
streaming. In its simplest form, CR have a Larmor radius 
much greater than the wavelength of spiral perturbations 
in a zeroth order uniform magnetic field. Because of their 
large Larmor radius the streaming CR, carrying an elec- 
tric current density ] CR , are essentially undeflected by the 
perturbed field but the jcr x B force acts towards the cen- 
tre of the spiral. A corresponding reactive force acts on the 
background plasma to expand the spiral. This stretches and 
increases the magnitude of the perturbed magnetic field, 
thereby increasing the jcr X B force in a positive feed- 
back loop that drives the instability. NRH appears to be 
the most rapidly growing instability driven by CR stream- 
ing on a relevant scalelength. Other instabilities that have 
attracted significant attention are the resonant Alfven in- 
stability (Kulsrud & Pearce 1969, Wentzel 1974) that grows 
with spatial wavelengths spatially resonant with the CR Lar- 
mor radius and the Weibel instability (Weibel 1959) that 
grows quickly on the spatial scale of an electron or pro- 
ton collisionless skin depth, c/ui pe = 5.3(ii e /cm~ 3 ) _1 ' 2 km, 
c/ujpi = \J m p /m e c/ujpe. The Weibel instability is important 
for interactions engaging thermal or mildly suprathermal 
electrons and ions, and may be effective for CR acceleration 
to low energies, but it grows on a scale too small to scat- 
ter PeV ions which have a Larmor radius of ~ 10 9 c/ui pe . 
The Alfven instability grows on the desired spatial scale 
but grows less quickly than the NRH instability in the SNR 
shock environment. Instabilities that grow on scales larger 
than the CR Larmor radius (Bykov et al 2011, Drury & 



Falle 1986, Drury & Downes 2012, Malkov & Diamond 2009, 
Rogachevskii et al 2012, Schure & Bell 2011) tend to have 
longer growth times but turbulent amplification may in- 
crease the growth rate (Bykov et al 2011). An extended dis- 
cussion of instabilities driven by cosmic ray streaming can 
be found in Schure et al (2012). 

From the Lagage & Cesarsky (1983a,b) limit as de- 
scribed in section 1 it is clear that proton acceleration to 
a few PeV is only possible if the magnetic field is strongly 
amplified above its characteristic interstellar value of a few 
pG. The following argument places an upper limit on the 
maximum energy of a CR that can be strongly scattered by 
a magnetic field growing on the scale of a CR Larmor radius. 
The background thermal plasma is highly magnetised on the 
scale of a PeV CR Larmor radius. Consequently the mag- 
netic field is 'frozen in' to the thermal plasma. Magnetic field 
amplification occurs as the plasma moves and stretches field 
lines. Assuming the perturbed magnetic field does not far 
exceed the initial field, the jcr x B force displaces a plasma 
element a maximum distance Smax ~ (jcRBt 2 )/2p in time t 
where B is the initial seed field. For the stretched magnetic 
field to strongly scatter CR as required for Bohm diffusion 
the displacement must grow to the order of a CR Larmor 
radius, that is s max ~ T/cB where T is the CR energy in 
cV. The CR current density jcr carries an energy flux jcrT 
in the upstream plasma rest frame which cannot far exceed 
the energy flux pu^/2 carried by the upstream plasma into 
the shock, where p is the upstream mass density. We de- 
fine a CR acceleration efficiency r\ such that jcrT = rjpu 3 a . 
We then have two equations: s ma x = (jcRBt 2 )/2p ~ T/cB 
and jcR = T)pUs/T. When combined, they yield a CR en- 
ergy T ~ (rjcul) 1 ^ 2 Bt. This expression for T is equivalent 
to T ~ 1.5(^) 1/2 B MG iioooPeV where U7 is the shock ve- 
locity in units of 10,000km s , B^g is the seed magnetic 
field in pG from which amplification begins, and tiooo is the 
time in 1000's of years. Characteristically for young SNR, 
ur — 0.5, tiooo = 0.4, r\ = 0.03 (see Appendix A) and 
£>mG = 3 where the magnetic field of a few pG represents 
the seed field from which the instability grows. With these 
estimates, T ~ O.lPeV, which is an order of magnitude less 
than the energy of the knee. This estimate is independent 
of any detailed instability theory except for the assumption 
that magnetic field amplification takes place through plasma 
motions generated by the }cr X B force acting on the scale 
of a CR Larmor radius. It highlights the difficulty of ampli- 
fying magnetic field by orders of magnitude on the scale of 
the Larmor radius of a PeV proton and the need for an in- 
stability that can provide non-linear growth of the magnetic 
field. 

The NRH instability offers a way round this diffi- 
culty by initially growing very rapidly on a scale much less 
than the CR Larmor radius. Since the small-scale mag- 
netic field grows unstably the jcfl X B force grows expo- 
nentially in time. The scale-size increases to the CR Lar- 
mor radius during non-linear growth. The NRH growth 
rate is 7 = (kjcRBo/p) 1 ^ 2 where Bo is the zeroth or- 
der magnetic field and k is the wavenumber on which 
the instability grows. The maximum NRH growth rate is 
7™i = 0.5 jcR^o/p) 1 ' ' 2 which occurs at a wavenumber 
k m ax = 0.5p jcR/B . "f max is (k max r g ) 1/2 times larger than 
the NRH growth rate for kr g = 1, thus allowing the mag- 
netic field and the jcr X B force to increase rapidly. The 
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NRH and Alfven growth rates are similar at kr g = 1 (see 
Appendix E and Bell 2004). Using the above nomenclature, 
kmaxTg = r)u 3 s /cv 2 A = 2 x 10 4 77o.o3M7W e B~Q where va is the 
Alfven speed, r\ = 0.03r;o.o3 and n e is the electron density 
in cm -3 . For a discussion of the value of r) see Appendix 
A and Bell (2004). The growth time of the fastest growing 



mode is then 7 maa; = 50% 



03 



in, 



-!/2„ 



Tp e y years where 



Tpev is the energy in PeV of the CR driving the instability. 
"/max — 400 years for ??o.o3 = 1, n e = 1, 117 = 0.5 which 
implies that even the NRH instability struggles to amplify 
the magnetic field sufficiently to accelerate CR to PeV en- 
ergies in the historical SNR. For growth by many e-folding 
the NRH instability must be driven by CR with energies less 
than IPeV. 

As remarked above, the fastest growing mode grows on 
a scale k m \ x whereas efficient CR scattering requires fluctu- 
ations in the magnetic field on a scale r g . The above analysis 
showed that k m a X r g — 2 x 10 4 rio.o3U7n e B~Q. The instabil- 
ity initially grows on a scale too small to effectively scatter 
PeV CR. Two factors save the situation. Firstly, as the mag- 
netic field grows from a few fiG to a few hundred ^tG the 
CR Larmor radius decreases by the corresponding factor of 
~ 100. Secondly, the characteristic scalelength of the struc- 
tured magnetic field increases during non-linear growth of 
the instability. Figure 1 is drawn from figures 3 and 4 of 
Bell (2004). The graph of energy densities shows how the 
magnitude of the magnetic field grows to many times its 
seed value in a time ~ Wj^ax- The characteristic scale- 
length grows during this time as shown by the grey-scale 
images. By the time t = 10 the scalelength has grown to the 
size of the computational box from its initially much smaller 
scale. Both k and r g decrease by a large factor by the time 
Wjmax- The simulation reproduced in figure 1 had a fixed 
CR current so it was impossible to demonstrate strong CR 
scattering by the magnetic field after t = 107^^, but the 
figure provides strong evidence for rapid initial growth of the 
fastest growing modes into large scale magnetic structures 
able to scatter CR on the scale of a Larmor radius. 



3 CR ESCAPE UPSTREAM OF THE SHOCK 

The previous section does not provide an estimate of the 
magnitude of the amplified magnetic field, but it does de- 
fine conditions under which strong magnetic field amplifica- 
tion occurs. CR acceleration to PeV energies requires strong 
magnetic field amplification which fixes the number of insta- 
bility e-foldings to the range 5-10. This in turn fixes the CR 
current. From figure 1 we take the condition for strong mag- 
netic field amplification in a particular volume of upstream 
plasma to be that J "/maxdt ~ 5 in the time before the shock 
overtakes it. Since ^max = Q-^3CRyPo/p, the condition for 
field amplification is 



ICR 



= j jcadt = 10^ p/p 



(1) 



where Qcr is the total electric charge of CR passing through 
unit surface area upstream of the shock before the shock 
arrives. (It makes no difference whether the CR passing 
through the plasma have high or low energies. It is only 
the number of CR escaping upstream that counts.) Since 
only the highest energy CR escape, the areal charge Qcr 



is carried by the highest energy CR being accelerated. The 
Larmor radius of a PeV proton in the interstellar medium 
is comparable with the radius of the historical SNR and 
its scattering mean free path is very much greater than the 
SNR radius. These escaping CR pass into the interstellar 
medium with low probability of further encounter with the 
SNR shock. Strongly scattered CR at lower energies are con- 
fined by the shock and are advected away downstream after 
acceleration. 

Suppose that the CR distribution follows a p~ 4 power 
law up to a momentum p m ax, then in steady state the elec- 
trical current of CR escaping in a small band of energies 
above the energy eT max 



jCR = eu s np max fo(p m ax) 



(2) 



where fo is the isotropic part of the CR distribution in mo- 
mentum space at the shock. Equation (2) is derived from 
equation (11a) below by integrating in space across the shock 
and deriving the rate at which CR reach the momentum 
Pmax- The derivation can be found in Appendix B. In com- 
parison the CR pressure at the shock, where the CR distri- 
bution extends down to p = mc and m is the proton mass 
is 

4tt 



Pcr = -ircp rnax fo(pmax)\n 

3 \ mc 



ax \ 
,C ) 



giving 



j C R = 0.05 



Us Pcr 

Tmax 



(3) 



(4) 



where we have assumed \n(p ma x/mc) = 14 for a CR dis- 
tribution extending from lGeV to IPeV. The condition for 
magnetic field amplification that f jcRdt = 10-y/ p/po can 
be rearranged to give a value for T max : 



T 



= 0.005- 



which is equivalent to 



o 1/23, PcRt, -, r 
8 V u 7 tiooo 5- PeV 



(5) 



(6) 



For characteristic values (117 = 0.5, n e = 1, t — 0.4, 
Pcr/puI = 0.3) T ma x ~ lOOTeV. As expected from the 
above discussion, this falls short of the few PeV required 
to explain Galactic CR. This will be discussed further in 
section 8. 

Zirakashvili et al (2008a,b) developed a related analytic 
model of the excitation of the NRH instability and CR con- 
finement upstream of the shock. They derived an estimate 
for the maximum CR energy with a similar dependency on 
nl^ 2 u 3 t in the numerator of their equation (19) but with 
an additional denominator that depends on the magnitude 
of the amplified magnetic field. Their analysis operates at a 
more detailed level than ours since they consider the growth 
in amplitude of the amplified magnetic field over a range 
of scales and small angle CR scattering by small scale field. 
However, the dominant physics is similar in their analysis 
and ours, and their estimate of the CR energy (equation 
(19)) is similar to that in our equation (6). 

Note that the magnetic field does not enter into the es- 
timation of the maximum energy CR in equation (6). The 
magnetic field is assumed to grow to whatever magnitude 
and spatial scale required to confine CR at energies less 
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Figure 1. Plots of the magnitude of a magnetic field driven by CR streaming into the page (2D slices of a 3D simulation). The field 
grows from noise on a small spatial scale at t = 0. By t = 10 the spatial scale has grown to the size of the computational box. The graph 
shows how the mean magnetic, kinetic and thermal energy densities increase with time. This figure is composed from figures 3 and 4 of 
Bell (2004) where further details can be found. The units of time are such that -f m ax = 1.26. 



than T m ax- This assumption is justified on the basis that 
the magnitude of the field increases by a substantial numer- 
ical factor in time 7~ nl even in the non-linear regime and 
that this is accompanied by rapid growth in characteristic 
spatial scale. For example, the magnetic field and its spatial 
scale are very much larger in figure 1 when J ^maxdt m 8 
than when J *y m axdt ~ 5 even though the time is different 
by a factor of only 1.6. Allowing the field to grow for just 
a little longer gives a much larger magnetic field and spa- 
tial scale. Consequently the primary parameter is not the 
magnitude of the magnetic field but the charge carried by 
escaping CR. 

It is clear that magnetic field amplification is only possi- 
ble if a population of high energy CR escape upstream of the 
shock. A magnetic field capable of stopping their escape is 
only produced after the areal charge Qcr has escaped. Thus 
CR escape upstream of the shock is an essential aspect of 
diffusive shock acceleration when magnetic field amplifica- 
tion is operational. If a charge Qcr has not passed through 
a particular point upstream of the shock then CR cannot be 
confined at that distance from the shock. 

Previous authors have also concluded that CR escape 
upstream is inevitable, but generally for different reasons. 
One possibility is that CR escape through filamentary cav- 
ities in the magnetic field (Reville & Bell 2012, 2013) but 
usually it is assumed that CR escape at a free escape bound- 
ary in position or momentum. Free escape boundaries were 
introduced at an early stage in the development of shock 
acceleration theory (Ellison et al 1981). It was appreciated 
that steady state models were impossible without CR escape 
(Ellison & Eichler, 1984, Berezhko & Ellison 1999, Malkov 
et al 2000). The standard T~ 2 test particle energy spectrum 
at strong shocks diverges if integrated to infinite energy. The 
divergence is even stronger if non-linear CR feedback onto 



the shock structure is included since the relativistic CR pres- 
sure increases the density jump at the shock resulting in 
a spectrum flatter than T~ 2 at high CR energy. However, 
the conditions for a steady state solution do not present a 
compelling argument for a free escape boundary since time- 
dependent solutions offer an acceptable option. Indeed the 
Lagage & Cesarsky limit on CR energy is based on the as- 
sumption that the maximum CR energy increases with time. 
The case for CR escape was strengthened by the recognition 
that magnetic field amplification is essential since there will 
always be a distance upstream at which amplification is in- 
operative and CR must escape. Similarly, there must always 
be an upper limit to the spatial scale to which field is ampli- 
fied and a corresponding limit, through the Larmor radius, 
on the momentum to which CR can be accelerated. This 
led to the imposition of a free escape boundary either at an 
imposed distance upstream of the shock (eg Caprioli et al 
2010b, Ohira et al 2010) or at an imposed CR momentum (eg 
Ellison & Bykov 2011). Whether the free escape boundary 
is imposed in momentum or configuration space, the mo- 
mentum at which CR escape depends on the magnetic field. 
Some authors (eg Ptuskin et al 2010) assume that the field is 
amplified until the magnetic energy density reaches a given 
fraction of pu 2 B . The fractional magnetic energy density can 
be chosen to match observation (Ptuskin et al 2010, Volk et 
al 2005). Alternatively the field can be chosen to match the 
saturation value estimated by Bell (2004). Others (Caprioli 
et al 2010a, Drury 2011) choose a mathematical form for the 
amplified magnetic field that allows for multiple possibilites. 
Ptsukin & Zirakashvili (2003) and Caprioli et al (2009a,b, 
2010b) include magnetic field generation due to the resonant 
instability as it develops ahead of the shock. Vladimirov et 
al (2006) include magnetic field generation in response to 
gradients in the CR pressure. Ellison et al (2012) choose an 
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amplified magnetic field suitable for SNR expansion into a 
circumstellar wind. 

In contrast to the above, CR escape in our model is 
determined by the escaping CR electric charge rather than 
the magnetic field generated in the upstream plasma, al- 
though magnetic field amplification is implicitly required by 
our model. Reville et al (2009) also used the escaping flux to 
calculate growth rates in an effort to motivate a realistic free- 
escape boundary location in their steady-state non-linear 
model. Their discussion on self-regulation of CR precursors 
however did not extend to the maximum CR energy. 



4 A NUMERICAL MODEL 

We now set out to test the above conclusions as far as 
we are able with a numerical model that includes the 
self-consistent interaction of CR modelled kinetically with 
a background plasma modelled magnetohydrodynamically. 
Standard MHD equations describe the background plasma 
except that a —]cr x B force is added to the momentum 
equation: 

P^ = - Vp -— Bx(VxB)-j C8 xB (7) 
at no 

as described in Lucek & Bell (2000) and Bell (2004). The CR 
distribution function /(r, p, t) at position r and momentum 
p is defined in the local fluid rest frame and evolves accord- 
ing to the Vlasov-Fokker-Planck (VFP) equation 

f t =-vM +Pi ^K- emeViB V +C (f) (8) 
at an dri dpj dpk 

where C(f) is an optional collision term included to repre- 
sent scattering by magnetic fluctuations on a small scale. 
The electric field is zero in the local fluid rest frame. 
Quadratic terms in the local fluid velocity u are neglected 
on the assumption that u <C c. The CR current jcr, re- 
quired for insertion into the MHD momentum equation, is 
calculated by integration over / in momentum space. The 
magnetic field in the CR kinetic equation is taken from the 
MHD calculation. 

As in Lucek & Bell (2000) and Bell (2004) three spa- 
tial dimensions are needed to represent the turbulence ad- 
equately. Since our aim is to investigate the mutual inter- 
actions of magnetic field amplification, CR acceleration and 
CR escape upstream of the shock we need to model the 
complete system including the shock and the complete CR 
precursor. Because we model the detailed interaction be- 
tween the CR and the distorted magnetic field we need to 
resolve the CR Larmor radius in configuration space and the 
rotation of CR trajectories in momentum space. As a con- 
sequence the numerical model should be 6-dimensional in 
momentum-configuration space with spatial scales extend- 
ing from the CR Larmor radius of the lowest energy CR 
to the precursor scaleheight of the highest energy CR. This 
would be impossible without extraordinary computational 
resources so our strategy is to design a computational model 
that includes all the important processes at a minimal level. 
We retain the three dimensions in configuration space but 
limit the range of CR momentum to a factor of 10 so that we 
do not have to resolve the Larmor radius of very low energy 
CR. We choose a shock velocity u 3 = c/5 to keep the ratio 



of the CR to the MHD timescale to a minimum while stay- 
ing close to the range of conceivable SNR expansion speeds. 
Our greatest approximations are made in the momentum 
space representation of the CR distribution function since 
this is the aspect of the calculation in which the number of 
dimensions can be reduced. 

The Vlasov-Fokker-Planck (VFP) equation (equation 
8) is important in the physics of laser-produced plasmas 
where it is solved in finite difference form to model electron 
transport. The successful use of VFP simulation to model 
electron transport in laser-produced plasmas stretches back 
more than 30 years (Bell et al, 1981) so it is natural to ap- 
ply the techniques to CR which obey the same equation. The 
distribution function /(r, p, t) of charged particles (cosmic 
rays or energetic electrons) is usually represented in spher- 
ical co-ordinates (p,8,(j>) in momentum space. A common 
representation of the distribution function is as a sum of 
spherical harmonics: 

fir, P ,t)=j2 /" i ( r - p> *) p « |m| ( cos e y m * 

I .m 

i = 0,oo m=-l,l /r m = (/D* (9) 

where /™(r,p, t) is the coefficient for the (l,m) spherical 
harmonic. / ; m (r, p, t) is function of time, position and mag- 
nitude of momentum p. The spherical harmonics decribe the 
angular structure on shells of constant magnitude of momen- 
tum. Reviews of the VFP technique, or papers containing a 
significant review element, are Bell et al (2006), Tzoufras 
et al (2011) and Thomas et al (2012). For information on 
the application of VFP techniques to CR acceleration the 
reader is referred to Reville & Bell (2013) which includes an 
appendix setting out the full VFP equations for a spheri- 
cal harmonic expansion. Bell, Schure & Reville (2012) also 
apply VFP techniques to CR acceleration. The use of an ex- 
pansion in tensors (used here) as an alternative to spherical 
harmonics is discussed by Schure & Bell (2012). 

VFP simulation was used by Bell et al (2011) for the cal- 
culation of CR acceleration by oblique shocks. They found 
that an expansion to 15th harmonic could be needed for 
oblique shocks because of the abrupt change in magnetic 
field direction at the shock, but that only a few harmon- 
ics are needed for quasi-parallel shocks. Here we reduce 
the computational size of the problem by modelling parallel 
shocks in which the zeroth order magnetic field is parallel 
to the shock normal. There are good reasons to suppose 
that the first few terms in the expansion capture the es- 
sential physics. Firstly, as shown by Bell (2004) the NRH 
instability is driven by the CR current density ]cr- Higher 
order anisotropics do not directly contribute to the insta- 
bility. Secondly, the CR precursor scaleheight is c/u s times 
the CR scattering mean free path in standard DSA theory. 
DSA theory is based on the diffusive approximation in which 
only the first order anisotropy is needed and only the first 
two terms (isotropy and drift anisotropy) in the harmonic 
expansion are retained. In diffusion theory the higher or- 
der anisotropies are damped by scattering. Hence it might 
be supposed that only the first two terms are needed and 
an adequate representation of the CR distribution func- 
tion might be /(r,p,t) = /o(r,p,t) + /°(r,p,t) cos0 + 
U{fi (r, p, t) sin 6>e^}. This '/o+/i ' expansion allows free CR 
propagation along magnetic field lines but restricts trans- 
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port across the magnetic field because the direction of the 
anisotropic drift term is rotated by the field. However, this 
expansion omits an essential feature of CR transport. It does 
not allow CR to gyrate as they travel along a magnetic field 
line. The fo + /i expansion allows CR to propagate along 
field lines and separately it allows CR to gyrate around field 
lines, but it does not allow CR to do both at the same time. 
Spiral trajectories require the inclusion of the off-diagonal 
components of the stress tensor fi. Without the stress ten- 
sor, CR cannot resonantly interact in space with magnetic 
perturbations on the scale of a Larmor radius. Clearly this 
would rule out the Alfven instability, and more surprising 
it also rules out the NRH instability. The role of the stress 
tensor is discussed by Schure et al (2011) in which the linear 
NRH dispersion relation is derived with the perturbed CR 
distribution expressed as a tensor expansion. We proceed 
on the basis that the CR distribution function may be ade- 
quately represented by an isotropic part (0th order in u s /c) 
plus a drift component (1st order in u s /c) plus a term rep- 
resenting the stress tensor (2nd order in u s /c). This second 
order expansion is more easily represented in the equivalent 
tensor notation instead of spherical harmonics: 

/(r, p, t) = f (r, p, t) + fi (r, p, t) * + f tj (r , p, t) ^ (10) 



V 



P 



where the trace of fij is zero because it is already accounted 
for in the isotropic term fo. Following Johnston (1960) the 
reduced Vlasov equation for the CR distribution function is 
then: 

dfo d(foUj) = cdfj duj 1 d(p 3 f ) _ , 
dt dri 3 dti dri 3p 2 dp 



dfi d(fiUj) 
dt dri 



dfo 2c dfij 



dri 



dfg 
dt 



+ 



d{fijU k ) 



dr k 



5 dr 

dfi 
dri 



ijh^fk (Hb) 



+ 



dfi 

dr t 



+ ^ 



' dr k 



ceB k 



(tkilflj + tkjlfli) 



(lie) 



The notation is simplified by introducing the vector gi 
related to shear and vorticity (V x fi) in CR motion to rep- 
resent the off-diagonal components of the stress tensor. The 
on-diagonal components of /y (i — j) are approximately 
accounted for by multiplying the —cdfo/dri term in equa- 
tion (lib) by 9/5 to allow for the stress tensor contribution 
to compressional waves and to allow freely streaming CR 
to propagate at y/iJEc instead of \JT~/3c, as shown in Ap- 
pendix D. The equations to be solved (their derivation is 
given in Appendix C) are then 



dfo 
dt 



+ 



d(foUi 



dr t 



cdfj | dm 1 d(p s f ) 
3 dri drt 3p 2 dp 



dh 
dt 



+ 



drj 

dgi 
dt 



9 dfo 1 dg k ceBj 

"7 c a r cei i k '" 

5 dri 5 



dr 



<Hjk - 



+ 



djgiUj 

dri 



dfh 
dri 



VBg l 



fh 



(12) 



Equivalently, expressed in vector notation, 

dfo c V.u gV/q) 

— + V .(u/o) - --V.fr + —-^- 



-^v/o-fvxgi-nxfi 

o 



ggl 

dt 



+ V.(ugi) = cV X fl - VB%\ 



(13) 



where ft = ecH/p is the vector CR Larmor frequency and 
we take vb = ecB/p. The presence of the curls of fi and 
gi facilitates the propagation of transverse modes in fi and 
gi as needed for helical motion along magnetic field lines or 
CR propagation at an angle to the wavevector k. 

The termination of the harmonic expansion at the stress 
tensor makes the computation tractable with available re- 
sources. Appendices D and E show that the truncated ex- 
pansion provides an adequate representation of the essential 
physics of CR propagation and CR-driven instability. 



where quadratic order terms in the velocity u have been ne- 
glected and we have omitted terms involving df /dp times a 
gradient of u apart from the term in equation (11a) for the 
evolution of fo. This amounts to the neglect of second or- 
der Fermi acceleration and acceleration resulting from shear 
motions in the background hydrodynamics. In our problem, 
these processes are small in comparison with acceleration at 
the shock. 

We have chosen to terminate the set of equations at 
equation (11c) by setting fij k to zero. Additionally, we soften 
the termination of the tensor expansion by replacing the 
magnetic rotation in equation (11c) by a damping term with 
a damping rate equal to the magnetic gyration frequency. 
The logic of this approximation is that random rotation of 
the stress tensor anisotropy leads to its damping. A similar 
assumption underlies Bohm diffusion which replaces rota- 
tion of fi by a damping rate, thereby terminating the tensor 
expansion at fi whereas we terminate it at /y . Our approx- 
imation makes intuitive sense, and we support it in Appen- 
dices D and E by examining its effect on propagating modes 
and on the NRH instability. 



5 THE SIMULATION 

A full 3D simulation with realistic parameters is not possi- 
ble because of the large ratio of the largest distances (the 
CR precursor scaleheight and the CR free escape distance) 
to the smallest distance (the shortest wavelength on which 
the NRH instability grows). The correspondingly large ra- 
tio of the largest timescale (the SNR expansion time) to 
the shortest timescale (the shortest NRH growth time) simi- 
larly makes substantial demands on computer resources. The 
computational constraints are discussed in Appendix F. 

We artificially increase the magnetic field and stretch 
other parameters in a favourable direction, choosing: 

n e = Clem -3 u s — 60,000km s" 1 



Bo = 47^iG v A = 3 X 10 5 m s T ln 



lOOTeV (14) 



where Ti n j ect is the energy at which CR are injected at the 
shock. The zeroth order magnetic field is aligned along the 
shock normal. The instability is seeded by imposing random 
fluctuations on the magnetic field with a mean magnitude of 
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9/j.G. CR are injected at the shock into the lowest momentum 
bin (width Ap) according to the rule 

2 ® fo 2 

47rp Ap c P~q£~ ~ —constant x V.u min(po , u , U) (15) 

where po is the density upstream of the shock and u is the 
local background fluid velocity relative to the initially sta- 
tionary upstream plasma. U is the local thermal energy den- 
sity. This prescription is designed to inject a suitable energy 
density of CR at the shock, dependent on the choice of the 
constant, whilst avoiding a negative value of U due to exces- 
sive transfer of energy to CR from cold plasma at the foot 
of the shock. The resulting CR current density ahead of the 
shock is displayed in panel (e) of figures 2 and 3. The peak 
CR current density in the population freely escaping ahead 
of the shock in figure 2 is 

j C R = 1.1 x lCT 14 Amp m" 2 (16) 

This current density corresponds to r\ — 0.026, giving 

7~L = 2.3 x 10 6 sec c 7m L = 7 x 10 14 m (17) 

kmlx = 7 x 10 n m r g = 7 x 10 13 m 

r g k max = 100 M A = 200 (18) 

We use a spatial grid with Aa; = Ay = Az = 1.4 x 10 12 m. 
Ten cells in momentum cover an energy range from lOOTeV 
to IPeV with logarithmic spacing. There are 32 cells in each 
of x and y with periodic boundary conditions. In contrast 
5676 cells are used in z with reflective boundary conditions 
for the MHD part of the code and for CR at the right 
hand boundary. CR reaching the left hand boundary are 
disposed of on the assumption that they escape freely. Cor- 
respondingly, the computational box extends 7.9 x 10 15 m by 
4.5 x 10 13 m by 4.5 x 10 13 m. The box-size in x and y is com- 
parable with the initial CR Larmor radius, but the Larmor 
radius contracts significantly as the magnetic field is ampli- 
fied. The cell-size is 7r _1 times the wavelength of the fastest 
growing mode, which is barely sufficient to represent the 
initial growth of the instability, but the characteristic scale- 
length of the instability increases rapidly as its amplitude 
grows. These parameters are only marginally sufficient to 
represent the physics, but for example a halving of the com- 
putational cell-size would increase the computational cost 
by a factor of 16. The simulation in figure 3 was run for 
6.2 x 10 7 sec = 2yr using 128 processors for 75 hours on 
the SCARF-LEXICON cluster at the UK Rutherford Ap- 
pleton Laboratory. The MHD part of the code is the same 
as that used by Lucek & Bell (2000) and Bell (2004). The 
VFP part of the code that models the CR uses a second or- 
der Runge-Kutta scheme in configuration space, a donor-cell 
scheme in magnitude of momentum, and the Boris algorithm 
for rotation in momentum due to magnetic field as used in 
particle-in-cell codes (Birdsall & Langdon 1985). The VFP 
equations are formulated in tensor notation rather than in 
spherical harmonics, but the two formulations are similar for 
the low order expansion used here, and their numerical so- 
lution is informed by experience with the KALOS spherical 
harmonic code (Bell et al, 2006) where Runge-Kutta advec- 
tion is found to be robust, sufficiently accurate, and a good 
fit to the form of the equations. Numerical diffusion due 
to limited spatial resolution is ameliorated by designing the 



simulation to initialise the upstream background plasma at 
rest relative to the computational grid. The more usual ap- 
proach of initiating the simulation by setting the background 
plasma in motion towards a reflective boundary would exac- 
erbate the effect of numerical diffusion on small structures 
during advection. Instead, a dense piston is initialised mov- 
ing leftward from the right boundary, pushing a shock before 
it into the stationary background plasma. 



6 SIMULATION RESULTS IN 3D 

Results from the 3D simulation are presented in figures 2 
and 3 at t = 4.1 x 10 7 sec and t = 6.2 x 10 7 sec respectively. 
Note that the horizontal spatial scales, in the direction of 
the shock normal, are artificially compressed by a factor 
of 24. The actual aspect ratio of the computational box is 
177:1. The plot of the background plasma density in panel 
(a) shows the position of the dense plasma piston propagat- 
ing leftwards and pushing the shock ahead of it. The high 
pressure region in panel (b) is due to plasma heating at the 
shock. The CR pressure is plotted in panel (c). The maxi- 
mum CR pressure occurs at the shock. Panel (c) of figure 2 
shows the CR population dividing into a population freely 
propagating ahead of the shock and a population confined 
by magnetic field (panel (d)) near the shock. CR confined 
near the shock continue to be accelerated. The banded struc- 
ture in the CR pressure ahead of the shock in panel (c) is 
exaggerated by the colour coding. It represents only a varia- 
tion of a few percent in the CR pressure and probably arises 
from a small oscillation in the injection rate at the shock. 
Also, the horizontal spatial compression of figure 2 distorts 
the aspect ratio of the bands. 

The escaping CR are an essential aspect of the process of 
CR confinement by the self-generated magnetic field. They 
excite the instability and carry the escaping areal charge 
Qcr = / jcndt identified in equation (1) as necessary for 
substantial field amplification. In figure 2 the escaping CR 
are only just reaching the left hand of the grid where the 
integral J jcndt is small and field amplification is negligible. 
In contrast, immediately in front of the shock in figure 2 
/ jcndt has reached the critical value needed for strong field 
amplification causing the onset of CR confinement. 

By the time of figure 3 most of the escaping popula- 
tion has passed through the free escape boundary at the left 
hand end of the grid. By this time a large magnetic field in 
the upstream plasma has switched off CR escape creating 
an expanse of low CR density between the confined and es- 
caping CR populations. Panel (d) shows that the magnetic 
field is amplified by an order of magnitude by the escaping 
CR over a large distance ahead of the shock. 

The separation of CR into escaping and confined pop- 
ulations matches expectations from the argument in section 
3. The CR charge per unit shock area in the escaping pop- 
ulation is ~ 1 x 10~ 7 Coulomb m -2 in figure 2 in good 
agreement with the estimate of 1.3 x 10~ 7 Coulomb m~ 2 
from equation (1). Panels (e) in figures 2 and 3 plot the CR 
current density ahead of the shock. 
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Figure 2. 2D slices of a 3D simulation with peak values in brackets: electron density (9.7cm -3 ), pressure (340nPa), CR pressure (2.3nPa), 
magnitude of the magnetic field perpendicular to the shock normal (610/jG), CR current parallel to the shock normal (1.1 X 10~ 14 Am~ 2 ). 
t = 4.12 X 10 7 sec (1.3 years). The dimensions of the computational box are 7.9 X 10 15 m by 4.5 X 10 13 m. The horizontal spatial scale, in 
the direction of the shock normal, is artificially compressed by a factor of 24. 
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Figure 3. 2D slices of a 3D simulation with peak values in brackets: electron density (6.9cm -3 ), pressure (320nPa), CR pressure (1.9nPa), 
magnitude of the magnetic field perpendicular to the shock normal (590/iG), CR current parallel to the shock normal (4.2 X 10 _15 Am — 2 ). 
t = 6.18 X 10 7 sec (2.0 years). The dimensions of the computational box are 7.9 X 10 15 m by 4.5 X 10 13 m. The horizontal spatial scale, in 
the direction of the shock normal, is artificially compressed by a factor of 24. 
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Figure 4. 2D simulations. Plots of f p 3 at CR energies 100, 167, 278 and 464TeV at times between 4.1 X 10 7 sec (top row) & 16 X 10 7 sec 
(bottom row), with computational box lengths L proportional to t between 7.9 X 10 15 m & 3.2 X 10 16 m chosen such that the CR travel 
the length of the box during the simulation time. The width of the computational box is 4.5 X 10 13 m in each plot. The vertical axes are 
in dimensionless units with peak values given against the axis. 



7 SIMULATION RESULTS IN 2D 

Figure 4 follows the calculations of figures 2 and 3 to later 
times with the same parameters but in 2D instead of 3D to 
reduce the demand on computer resources. The top row of 
plots in figure 4 is the 2D equivalent of the 3D results in 
figure 2 at t = 4.1 x 10 7 sec. The dimensions of the spatial 
grid are the same in figure 2 and the top row of figure 4. In 



subsequent rows of figure 4 the spatial dimensions normal 
to the shock are expanded by factors of 2, 3 and 4 such 
that CR travel the length of the grid in the respective times 
(t = 8.2 x 10 7 , 1.2 x 10 s & 1.6 x 10 8 sec). 

In 3D at t = 4.1 X 10 7 sec the dip in the CR density 
separating the escaping and confined CR has only recently 
formed as seen in figure 2(c). At the same time in 2D (figure 
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4 top row) the separation between the confined and escap- 
ing populations is relatively more developed. Otherwise the 
results are similar in 2D and 3D. The reason for the slightly 
earlier development of CR confinement in figure 4 is that 
magnetic structures can expand more freely in 2D. In 2D, 
magnetic structures expand in the ignored dimension with- 
out running into stationary or counterpropagating plasma. 
Nonlinear development in 2D is therefore less inhibited and 
the magnetic field grows more rapidly. The magnetic field in 
2D reaches a steady value similar to that found in the 3D 
calculation. The eventual limit on the magnitude of the mag- 
netic field may be set by magnetic tension which operates 
equally in 2D and 3D. 

Comparison of the CR distributions at different times 
in figure 4 shows that at early times (t — 4.1 x 10 7 sec) only 
the low energy CR are confined. At t = 8.2 x 10 7 sec CR 
are confined at 100 and 167TeV with confinement beginning 
to occur at 278TeV and no evidence of the development of 
two separate populations in the few CR reaching 464TeV. 
By t — 1.2 x 10 8 sec a larger number of CR have reached 
464TeV with a local minimum in CR density fop 3 evident 
at all energies up to 464TeV in front of the shock. By t — 
1.2 x 10 8 sec CR at 100 and 167 TeV are strongly confined 
with a short precursor scaleheight ahead of the shock, which 
is unsurprising since the Larmor radius of a lOOTeV proton 
in a 700/^G magnetic field is 5 x 10 12 m. By t = 1.6 X 10 8 sec 
CR are confined at all energies up to 464TeV. 

At all four times in figure 4 the CR spectrum is much 
steeper than the steady-state test particle spectrum (/o oc 
p~ 4 ) because of CR loss into the downstream plasma. Mag- 
netic field amplification takes place ahead of the shock and 
the downstream field only becomes large when the shock 
overtakes the field amplified in the upstream. Hence at 
early times CR escape downstream. This causes a reduc- 
tion in the number of CR accelerated to high energy at the 
shock, thereby steepening the spectrum. Downstream con- 
finement at the shock improves at later times as shown by 
the downstream gradients in fo at t = 1.2 x 10 8 sec and 
t = 1.6 x 10 8 sec. However, unphysical numerical relaxation 
of the spatially compressed downstream magnetic field due 
to limited spatial resolution may be playing a part in allow- 
ing CR to escape downstream. 

There is slight evidence of pulsed acceleration as seen 
in the three separate peaks in the 278TeV CR density at 
t = 1.2 X 10 8 sec and t = 1.6 X 10 8 sec but overall the CR 
profiles develop in an orderly manner. 



R, CR are confined if 



8 APPLICATION TO SNR 

Our simulations support the model developed in section 3. 
According to our model CR are confined and accelerated if 
the electrical charge of CR escaping upstream of the shock 
reaches Qcr = 10^/ p/po Coulomb m~ 2 . We now apply this 
to spherical SNR shocks. The CR current density at a radius 
R is jcR = r]pu 3 r 2 / R 2 T due to CR accelerated to energy eT 
when the shock radius was r. Since only the highest energy 
CR escape upstream we assume that the CR reaching the 
radius R are monoenergetic with energy eT. T evolves as 
the shock expands. When the SNR shock reaches the radius 



L 



R ^"fr) r 2 dr = 10R 2 
T(r) 



P(R) 
Po 



(19) 



Differentiating this equation with respect to R and assum- 
ing a power law dependence of density on radius, p{R) = 
p (R/R )- m , 



Defining 117 



T <*>-6(4-mr 
Us/10, 000km s" 1 , R p . 



RJ-p 



(20) 



??/0.03, n e — p/2 x 10 _21 kg cm -3 (such that n e is approxi- 
mately the electron density in cm" 3 ) 
expansion into a uniform medium, 



_R/parsec, 770.03 
hat n e is app 
and taking m = for 



T = 230 r)Q.0'inJ 2 u 7 R vc TeV 



(21) 



A SNR with 117 — 0.6, n e = 1 and R pc = 1.7, representative 
of Cas A, would then accelerate CR to ~140TeV. A SNR 
with U7 = 0.5, n e = 0.1 and R pc = 10, representative of 
SN1006, would accelerate CR to ~180TeV. These energies 
are a factor of 10 lower than the energy of the knee. Abbasi 
et al (2012) place the knee at 4-5PeV, although data from 
other experiments indicate a lower energy and the turnover 
in the spectrum is not well defined as shown in their figure 
15. 

Our model places a considerable question mark over 
the ability of the well-known historical SNR to accelerate 
CR to the knee. Acceleration by SNR such as Cas A and 
SN1006 fails to reach the knee in our analysis because their 
expansion is already significantly decelerated. Zirakashvili 
& Ptuskin (2008a) also reach the conclusion (their Table 2) 
that the historical SNR do not accelerate CR beyond 100- 
200 TeV, but note that their parameter rj esc differs from our 
parameter 77 by a factor of two (rj e3C = 2rf). 

Initial expansion velocities of very young SNR can reach 
30,000 km s" 1 (Manchester et al 2002) or possibly higher 
for some types of SN (Chevalier & Fransson 2006). Accord- 
ing to equation (21) expansion into a density of 1 cm -3 at 
30, 000 km s _1 for 16 years to a shock radius of 0.5 parsec 
accelerates CR to ~ IPeV. Moreover, the energy processed 
through the shock is comparable to that for expansion to 1.5 
parsec at 6, 000 km s _1 , thereby contributing a comparable 
CR energy content to the Galactic energy budget. A com- 
plementary perspective on the same problem is obtained by 
expressing the maximum CR energy in terms of the mass 
M — AnpR 3 /3 swept up by the shock and the characteristic 
energy of the blast wave E = Mu 2 s /2. If Mq is the mass in 
units of a solar mass and En is the energy in units of 10 44 J, 
the maximum CR energy is 



T = 0.5770.03"- 



1/6 



E 4i M~ 2/3 PeV 



(22) 



which indicates that a typical SNR should be able to ac- 
celerate CR to ~ IPeV and that the maximum CR energy 
is greater for the same energy E given to a smaller mass 
M. The maximum CR energy is nearly independent of den- 
sity n e because a shock expands to a larger radius in a low 
density medium before deceleration (Hillas 2006). 

In the self-similar Sedov phase, E is constant, M is pro- 
portion to _R 3 and hence T oc R~ 2 . The maximum CR energy 
decreases with radius during the Sedov phase until mag- 
netic field amplification ceases, at which point our analysis 
in terms of the charge carried by escaping CR is inapplicable. 
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Equation (20) with m — 2 indicates that the highest 
CR energies might be achieved by SNR shocks expanding 
into a dense circumstellar medium previously ejected as a 
wind from the SN progenitor. pR 2 is independent of radius 
in a steady wind, in which case the CR energy depends only 
on u 2 and r\ according to equation (20) thus favouring CR 
acceleration in the very early stages of rapid expansion as 
previously suggested by Volk & Biermann (1988) and Bell 
& Lucek (2001). The maximum CR energy for a SNR ex- 
panding into a wind carrying M5 x 10 -5 solar masses per 
year and shed with a velocity of V4 x 10 km s _1 is 



T = 760r/ .o:iU 2 \ 1 



TeV 



(23) 



indicating that PeV energies are attainable. CR may be ac- 
celerated to energies beyond the knee if the initial shock 
velocity is ~ 30,000 km s _1 and the shock expands into 
a particularly dense wind or otherwise dense circumstellar 
medium. 

The estimates made in this section are subject to con- 
siderable numerical uncertainty. For example our estimates 
for the efficiency r\ ~ 0.03 or the escaping charge Qcr = 
10-y/ p/po could be uncertain by a factor of 2 or 3. How- 
ever, the accumulated error in our estimates would need to 
be a factor of 10 for the well-known historical SNR to ac- 
count for acceleration to the knee. Our arguments are not 
completely watertight but we tentatively conclude that ac- 
celeration to the knee takes place in younger relatively un- 
decelerated SNR. 



9 THE CR ENERGY SPECTRUM 

In this section we discuss the energy spectrum of escaping 
CR integrated over the lifetime of the SNR. The energy T 
of escaping CR changes as the SNR evolves in radius and 
expansion velocity. The integrated energy spectrum of CR 
escaping into the interstellar medium (ISM) need not be the 
same as the spectra of CR at the shock during acceleration or 
that of CR carried downstream into the interior of the SNR. 
Related analyses based on different models of CR escape can 
be found in Caprioli et al (2010a), Drury (2011), Ohira et 
al (2010), Ptuskin & Zirakashvili (2003) and Ptuskin et al 
(2010). 

We assume a power-law density gradient, p oc R~ m 
where m = for a uniform circumstellar medium and m — 2 
for a steady pre-supernova wind. We further assume that the 
shock velocity can be approximated as a power law u s oc R~ q 
over a sufficiently large part of the SNR's evolution. In this 
approximation u s oc t~ q/(1+q) and R oc t 1/{1+q) . In the Se- 
dov phase q = 3/2. From equation (20) T oc i? 1 " 2 "?-™/ 2 . 

In a uniform circumstellar medium (m = 0) the energy 
T of escaping CR decreases during expansion if q > 1/2, 
which is equivalent to u 3 decreasing more rapidly than t -1 ^ 3 . 
If a SNR expands into a steady pre-supernova wind (m = 2) 
the maximum CR energy always decreases with time pro- 
vided u 3 decreases with time (q > 0) as expected. 

Let E(T)dT be the total energy given to CR above 
an energy T. By definition of r\ the CR energy flux escaping 
upstream is r/pUg, so E(T)dT = f^ 1 r]4nr 2 pu 2 s dr and 

jn 

E = ~ ^nR 2 pu 2 s (24) 



where CR with energy T escape when the the shock radius 
is R. From equation (20) 



To 



l — 2q — m/2 



where To 



VMQ 2„ r— 

-— rUo-RoVPo 

5(4 — m) 



(25) 

and po, 1*0 and To are the values of p, u s and T at a reference 
radius R = Ro. Manipulating these equations gives 



CR 



4it V p u 2 R% (T_ 
(2q - I + m/2)eT 2 \T 



where a 



4g + 2 
4q + m - 2 



(26) 



and Ncr(T) = E/eT is the CR differential spectrum in en- 
ergy. For SNR expansion into a uniform medium the CR 
spectral index is a uni f orm — (2q + l)/(2g — 1), and for ex- 
pansion into a wind the index is a W i nc i = (2g+l)/2g. During 
the Sedov phase m = and q = 3/2, giving 



CtSe 



(27) 



although the analysis only applies while the magnetic field 
is being amplified by the NRH instability. A slightly less 
rapid decrease in shock velocity, u s oc t~ ' 57 (q = 4/3), 
would reproduce the CR spectral index (a m 2.2) inferred 
for CR at their source at energies less than IPeV (Gaisser 
et al 1998, Hillas 2005). The spectral index of CR escaping 
in the Sedov phase, asedov = 2, is the same as that for test 
particle acceleration at a strong shock. There is no obvious 
reason why this should be so. 

We emphasise that this discussion and the derived spec- 
tral index a or otsedov applies only to CR escaping upstream 
from the shock during SNR expansion. A further population 
of lower energy CR are carried into the centre of the SNR 
where they reside until they are released into the interstel- 
lar medium when the SNR slows, disintegrates and dissolves 
into its surroundings. These lower energy particles lose en- 
ergy adiabatically as the SNR expands but they can be ex- 
pected to contribute most of the Galactic population of low 
energy CR. 

In section 8 we suggested that acceleration beyond the 
knee might result from high velocity expansion into a dense 
circumstellar medium. The CR population beyond IPeV is 
an uncertain mix of protons and heavy ions. The overall 
spectral index of CR released into the Galaxy by SNR be- 
yond IPeV can be approximated as a w 2.7 with further 
spectral steepening occurring during propagation from the 
source to the Earth. This spectral index at source is pre- 



dicted by our analysis if q = 0.29 (u s oc t 



for expan- 



sion into a steady wind. It would correspond to a reduction 
in the shock velocity by a quite reasonable factor of 2.1 be- 
tween t — 10 years and t = 300 years. Expansion into a 
uniform medium would require q — 1.09 (u s oc t~ ' 52 ) equiv- 
alent to a reduction in shock velocity by a less reasonable, 
but not impossible, factor of 5.9 between t = 10 years and 
t — 300 years. This supports the contention that accelera- 
tion to 10 — lOOPeV may occur at a SNR shock expanding 
into a circumstellar wind. 
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10 MAGNETIC FIELD 

The Hillas parameter £u a BR (shock velocity times mag- 
netic field times spatial scale) provides an estimate of the 
energy T max to which CR can be accelerated under various 
circumstances (Hillas 1984). £ is a numerical factor of or- 
der unity, probably lying in the range between 1/8 and 3/8 
depending upon the CR diffusion coefficient (Lagage & Ce- 
sarsky 1983a,b, Bell 2012). T max = 4£ 8 M 7 B M G-RpcTeV where 
£ = Cs/8 (£s = 1 for £ = 1/8). Combining this with equa- 
tion (21) gives an estimate of the pre-shock magnetic field 
required to accelerate CR to T ma x'- 



the magnetic field given by equation (28) to the saturated 
magnetic field is 



B ~ 6077o.o3£g 1 n 1 J 2 u 7 pG . 



(28) 



The post-shock field is ~ 3 x larger due to compression at the 
shock. For our simulation parameters (n e = 0.1, ur = 6) the 
estimated post-shock maximum magnetic field is ~ 400/iG 
which is consistent with the field of ~ 600pG seen in panel 
(d) of figures 2 and 3. This confirms that the magnetic field 
in the simulation is amplified sufficiently to confine and ac- 
celerate CR and that diffusion in the CR precursor is ap- 
proximately Bohm, D ~ r g c. 

The above estimate of the magnetic field is derived on 
the basis that the magnetic field must be sufficient to con- 
tain and accelerate CR to the energy estimated in equation 
(21). Continued growth in the magnetic field would inhibit 
CR escape and remove the CR current that drives the NRH 
instability. Part of the energy generated by the NRH in- 
stability is stored in the kinetic energy of plasma motions. 
These motions might continue to stretch magnetic field lines 
and further increase the magnetic field after the CR cur- 
rent becomes inhibited. Further release of CR into the up- 
stream plasma would then be heavily restricted until the 
magnetic field relaxes to a lower level. This might result in 
oscillation about a marginal state defined by a balance be- 
tween magnetic field amplification and CR escape upstream. 
Weak evidence for periodic releases of CR into the upstream 
plasma can be found in the plot of the 278TeV CR density at 
1.2 x 10 s sec in figure 4 but on the whole the system appears 
to evolve without oscillation. 

In planar geometry escaping CR are in principle capa- 
ble of generating magnetic field at an unlimited distance 
ahead of the shock. In the spherical geometry of an expand- 
ing SNR the CR current decreases with distance ahead of 
the shock, jci? oc R~ 2 , so continuous CR escape is needed 
to amplify the magnetic field at a general radius R before 
the shock reaches that point. Hence the marginal balance 
between CR escape and magnetic field generation is more 
likely in spherical than planar geometry. 

The above discussion assumes that magnetic field 
growth and CR acceleration is determined by the growth 
rate of the NRH instability. However, it is possible that 
the instability might saturate and stop growing before it 
reaches that given by equation (28). Bell (2004, 2009) ar- 
gues that tension in the field lines limits amplification when 
V x B ~ pojcR for magnetic field structured on the scale 
of a CR Larmor radius. This implies a saturation magnetic 
energy density B 2 at /po ~ JcrT/c and predicts a saturated 
upstream magnetic field 



B s 



160 % 



1/2 1/2 3/2 



pG 



(29) 



B 

B sa t 



n Ai —\ 1/2 -1/2 



(30) 



2 2 



(31) 



with a further ~ 3x increase at the shock. The ratio of 



which implies that tension in the magnetic field does not 
stop the field growing to that given in equation (28). How- 
ever, it would not require magnetic field growth to overshoot 
that indicated by equation (28) by a large factor before the 
magnetic tension intervenes to halt growth. According to 
equation (28) the magnetic energy density is proportional to 
pu 2 , whereas the magnetic energy density determined by sat- 
uration is proportional to pu 3 s . Observations of SNR slightly 
favour a dependence on pu 3 s , but the difference is too close 
to be called (Vink 2006). 

From equation (20) (T = 0.05r)u 2 R(pp ) 1/2 for m = 
0) and the modified Hillas condition (T = u s BR/8) the 
required post-shock magnetic energy density is 

B^_ 

2p 

For magnetic field structured on the scale of the Larmor ra- 
dius of the highest energy CR we should assume r\ ~ 0.03 as 
the fraction of pu 2 given to the highest energy CR in which 
case the the post-shock magnetic energy density would be 
~ 0.1% of pu 2 , allowing for compression at the shock. Volk 
et al (2005) find observationally that the post-shock mag- 
netic energy density is typically ~ 3% of pu 2 in the histor- 
ical SNR. However, magnetic field will also be amplified on 
the scale of the Larmor radii of low energy as well as high 
energy CR. The difference between ~ 0.1% and ~ 3% may 
be explained by integration over the magnetic structures on 
scales varying by six orders of magnitude corresponding to 
the difference between the Larmor radii of GeV and PeV 
protons. If this is the case, most of the magnetic energy at 
the shock resides at scalelengths too short to accelerate CR 
to PeV energies. The magnetic field inferred from x-ray syn- 
chrotron observations of a SNR shock should not be inserted 
without adjustment into the Hillas parameter to estimate 
the maximum CR energy. 

Our analysis of CR escape leads to the result that the 
energy density of the magnetic field confining the highest 
energy CR is proportional to pu 2 . Previous analyses of CR 
escape often start from the assumption that the magnetic 
energy density is proportional to pu 2 and consequently they 
produce similar results for the spectrum of escaping CR. 
Comparable results for the CR spectrum produced in the 
Sedov phase can be found in Caprioli et al (2010a), Drury 
(2011), Ohira et al (2010), Ptuskin & Zirakashvili (2003) 
and Ptuskin et al (2010). For example a T~ 2 energy spec- 
trum for escaping CR in the Sedov phase has previously 
been derived by Berezhko & Krymskii (1988), Ptuskin & Zi- 
rakashvilii (2005), Caprioli et al (2010a) and Drury (2011). 



11 CR ENERGY INPUT TO THE GALAXY 

In our model only the highest energy CR escape upstream of 
the shock. At any given time more CR energy is carried away 
downstream into the SNR than escapes into the Galaxy. 
Efficient production of Galactic CR might therefore seem 
impossible. However, the CR energy carried into the interior 
of the remnant is subsequently recycled to drive the SNR 
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expansion and is available to accelerate further CR at a later 
stage. The overall efficiency of the production of Galactic CR 
is demonstrated by integrating over the lifetime of the SNR. 
Assuming a T~ 2 energy spectrum (a = 2, q = 3/2) equation 
(26) can be integrated to deduce a total CR energy input to 
the Galaxy: 

377 47^ 2 (T 2 \ 
E to tai = — —^—pou In I — I (32) 

between energies Ti and Ti (~ fGeV and ~ IPeV respec- 
tively). A small value of t; (~ 0.03) is balanced by the factor 
\u(T-ijT\) which results from recycling CR energy carried 
into the interior of the SNR. 

It is occasionally remarked that CR acceleration by very 
young SNR during the first few years cannot inject sufficient 
energy into the Galaxy to account for CR at energies beyond 
the knee because the SNR shock is small. However, if the CR 
spectrum connects smoothly across the knee and the spec- 
trum beyond the knee of escaping CR matches observation 
as shown to be possible in section 9 then it follows that the 
CR energy input is sufficient to match the Galactic energy 
budget. 



12 CONCLUSIONS 

The central message of this paper is that CR of a given 
energy escape freely ahead of a shock until magnetic field 
amplification takes place to inhibit propagation. The con- 
dition for propagation inhibition is that a sufficient num- 
ber of CR must escape upstream for the NRH instability 
to grow through ~ 5 — fO e-foldings at the growth rate of 
the fastest growing mode. Since the instability is driven by 
the CR current, the condition is that a CR electric charge 
Qcr ~ Wy/p/fio per unit area must escape through a spher- 
ical surface surrounding the SNR to amplify the magnetic 
field and inhibit CR escape through that surface. Since high 
energy CR carry less charge than low energy CR for a fixed 
CR energy flux the condition on Qcr determines the en- 
ergy of escaping CR. We find that the energy eT of escaping 
CR is proportional to r]Ru 2 a yfp as given by equation (21). 
The energy eT varies during the evolution of the SNR and 
determines the energy spectrum of CR injected into the in- 
terstellar medium by SNR. In our estimation, the historical 
SNR (Cas A, Tycho, Kepler, SN1006) are currently acceler- 
ating CR to ~ 100 - 200 TeV. Acceleration to the knee at 
a few PeV takes place in SNR at an earlier stage of evolu- 
tion when the shock velocity is ~ 10, 000 km s _1 or greater. 
This is an unsurprising conclusion since if the historical SNR 
were to accelerate CR to the knee we would be asking why 
even higher energy CR were not being produced by younger 
SNR. Observations by the planned Cherenkov Telescope Ar- 
ray (CTA) should be crucial in testing our conclusions (Hin- 
ton & Hofman 2010, Aharonian 2012). 

Acceleration beyond the knee may take place in very 
young SNR expanding at 20 — 30, 000 km s" 1 into a dense 
circumstellar pre-supernova wind. 

The spectral index of escaping CR is consistent with the 
measured Galactic CR spectrum at energies less than IPeV. 
Beyond the knee the proton spectral index is uncertain both 
theoretically and observationally. The theoretical prediction 



depends on the rate at which the SNR shock decelerates 
during its early expansion. 

The magnetic field can be estimated from the Hillas pa- 
rameter as the field needed to accelerate CR to the escape 
energy. The field is close to, but slightly less than, the satu- 
ration field determined by tension in the magnetic field. The 
predicted magnetic fields are consistent with those observed 
in SNR if allowance is made for the large range of scale- 
lengths, corresponding to the range of CR Larmor radii, in 
the structure of the magnetic field. 

The model is tested against numerical simulation. The 
MHD/VFP code is three-dimensional in space and one di- 
mensional in CR momentum with anisotropy included to 
second order. The computational parameters are pushed 
to their limit to allow solution of this multi-scale multi- 
dimensional problem but the results support the analytic 
model. In the simulation CR are seen to escape upstream 
with the electrical charge predicted by theory and magnetic 
field is strongly amplified to the predicted level. 
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Appendix A: The value of rj 

The CR electric current jcr drives the amplification 
of magnetic field through the NRH instability. Throughout 
this paper we express jcR as a fraction of the CR current 
needed to carry the characteristic energy flux pu^: jcn = 
rjpUg/T where T is the characteristic CR energy in eV. In 
this Appendix we briefly explain why we choose rj « 0.03 as 
our best estimate (see also Bell 2004). 

In the absence of CR acceleration the thermal energy 
density downstream of a strong shock is 9pu^/8 from the 
Rankine-Hugoniot relations. Assuming that a third of this 
energy is given to CR as required for efficient CR produc- 
tion by SNR, the downstream CR energy density is 3piig/8. 
From continuity across the shock the CR energy density im- 
mediately ahead of the shock is also 2>pu 2 a /8 and the CR 
energy flux relative to the upstream plasma is 3pUg/8. How- 
ever, only the highest energy CR escape upstream. Lower 
energy CR do not penetrate far upstream and they amplify 
magnetic field on too small a scale to engage the escaping 
CR. Hence the analysis in this paper depends on the cur- 
rent carried only by high energy CR. For a T~ 2 CR energy 
spectrum extending from T min a; lGeV to T max £s IPeV 
the energy is spread equally across each decade in energy 
with a fraction 1/ \n{T m ax/T min ) « 1/14 associated with 
any energy T. The energy flux carried by CR with energy T 
is then ~ 3pu 3 s /8ln(T max /T min ) ~ 0.03pu 3 s . The energy flux 
carried by CR streaming at velocity v with number density 
ncR and energy eT is ncRveT = jcrT. Consequently we 
assume that jcR = rjpu^/T where rj ~ 0.03. 

Appendix B: The derivation of equation 2 

Equation 2 for the electric current jcR carried by es- 
caping CR assumes (i) that CR escape upstream in a small 
range of momenta just above p ma x (ii) that any CR reach- 
ing a momentum p ma x or higher freely stream ahead of the 
shock at a velocity much greater than the shock velocity 
« s (iii) that escape is predominantly upstream. The electric 
current density is 

jcR = / ^-.ficp 2 dp (Bl) 

where fip r /p is the anisotropic part of the distribution fuc- 
ntion representing CR drift in the r direction. In a steady 
state on an acceleration timescale dfo/dt can be neglected 
from equation 11a. The second term d(foUi) / 'dri can also be 
neglected since the CR are assumed to be freely streaming 
at energies above cp ma x and advection at the fluid velocity is 
small. In one spatial dimension r equation 11a then reduces 
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to 



c dh , du 1 d(p 3 f ) _ Q 
3 9r dr 3p 2 dp 

for energies above cp ma x- Integrating in space across the 
shock and in momentum from p m ax upwards gives 



. 4ty 3 , . 

JCR = eAll—p fo{Pmax) 



(S3) 



where Aw is the change in velocity across the shock. Since 
An = 3« s /4 for a strong shock 



jCR = eirU s p f (Pmax) 



(B4) 



as in equation 2. 



Appendix C: The derivation of equations 12 and 13 

Here we show how equation 12 and its vector equivalent 
can be derived from equations 11. For simplicity we omit 
terms representing advection at the fluid velocity and the 
effect of the magnetic field. These can easily be inserted at 
the end of the derivation. The difficult part is the replace- 
ment of the stress tensor by the vector gi. The most 
transparent way of presenting the derivation is to write it 
out in terms of individual components in x, y and z di- 
rections. In the following we present the derivation of the 
equation for dfxjdt. The derivations of the equations for 
dfy/dt and df z /dt follow the same pattern. In component 
form, the relevant equations 11 are 

dfo _ c ( df x dfy dfz 
dt 3 \ dx dy dz 



These equations become equations 12 and 13 with the addi- 
tion of fluid advection, rotation of fi by the magnetic field, 
and damping of the stress tensor at a rate i/j as discussed 
above equation 12 in section 4. 

Note the analogy of fi and gi in the above equations 
with E and B in Maxwell's equations. In both cases they 
support transverse waves. 

Appendix D: The approximate CR equations: prop- 
agating modes 

Here we demonstrate that the equations set out in sec- 
tion 4 describe the essential CR propagation modes for mon- 
energetic CR. For simplicity we neglect the damping term 
(ys = 0). Propagating CR solutions in a stationary back- 
ground plasma can be found by setting u = 0, d/dri — > iki 
and d/dt — > —ica in equation (13): 

c 

9c c 
w/i = —kifo + -eijkkjgk — itijkQjfk 
o 



ojg m 



(Dl) 



The resulting dispersion relation is 

3c 2 fc 2 \ f ( , c 2 fc 2s 2 



: 



2fc 2 c 2 



(D2) 



dfx = _ c dfo _ 2c (df_ 
dt dx 5 \ dx 



dy 



dz 



dfxx = £ f 2 dfx dfy df z 
dt 3 \ dx dy dz 



dfxy c (dfx dfy 
dt 2\dy dx 



dfx 



dt 



- ( 2Jk _i_ d/z 
'2 V dz dx 



(CI) 



Eliminating the components of the stress tensor between 
these equations gives 

d 2 f x 9c d_dfo_ = l_ ( d 2 u + &% + d 2 fx 



(C2) 



dt 2 5 dt dx 5 \ dx 2 dy 2 dz 2 

_Sd_ (dfx + dfy + df 1 

5 dx \ dx dy dz 
which is the x component of the vector equation 



a 2 fi , 9c9(V/ ) 



dt 2 + 5 



-V x (V x fi) 



(C3) 



dt 5 

and the same derivation holds for the y and z components 
of the equation. This equation is second order in time dif- 
ferential. It can be separated into two first order equations: 

dfi 9c, 



dt 5 :V/o ~ 5 V X Sl 



dgi 
dt 



cV x fi 



(C4) 



where 9 is the angle between k and B. There are two inde- 
pendent modes in the absence of magnetic field (fi = 0). The 
mode propagating at ^/3/5c represents the motion of freely 
propagating CR with fi parallel to k. The mode propagating 
at yf\]Zc is the transverse mode representing the motions of 
CR with fi perpendicular to k. The transverse mode propa- 
gates more slowly because CR velocities are aligned prefer- 
entially away from the direction of propagation. 

In the presence of a magnetic field a longitudinal mode 
still propagates parallel to B (9 = 0) at ^/3/5c representing 
free propagation along field lines unaffected by the field. The 
transverse mode is relatively unaffected by the magnetic field 
at wavelengths shorter than a Larmor radius kc ^> fi. At 
wavelengths longer than the Larmor radius the transverse 
mode propagates more slowly as the transverse CR current 
rotates rapidly in the magnetic field. 

When mode propagation is directed across the magnetic 
field (sin^ = 1) the wave frequency is given by 



2 2 J 1 
u) — fi < - 



2k 2 ri 



^4 + 



25 



9 + ^ > (D3) 



lengths smaller than the Larmor radius, the frequency con- 
verges to those derived for the longitudinal and transverse 
waves in zero magnetic field as expected. At long wave- 
lengths (kr g <<iC 1) the frequency converges to the Larmor 
frequency, representing CR rotation in the magnetic field, 
without significant propagation. 

This analysis of the dispersion relation indicates that 
equations (13) for fo, fi and gi provide an adequate repre- 
sentation of CR propagation. 



© 2013 RAS, MNRAS 000,[2fj7] 



16 A.R. Bell, K.M. Schure, B. Reville and G. Giacinti 



non-resonant polarisation 
tensor expansion 



resonant polarisation 
tensor expansion 



102 , 4 



10° 



1° 2 kr„ 1° 4 



non-resonant polarisation 
Bell(2004) 




resonant polarisation 
Bell(2004) 



Figure 5. Dispersion relation for the resonant and non-resonant 
circular polarisations as derived from the tensor expansion (top 
row) compared with the dispersion relation derived by Bell (2004) 
for mono-energetic CR (bottom row). The growth rates in units 
of the CR Larmor frequency are given by full lines and the real 
frequencies by the dashed lines. Parameters relevant to the his- 
torical SNR are assumed. 



Appendix E: The approximate CR equations: the 
NRH instability 

We now investigate whether the approximate treatment 
of CR kinetics in section 4 is adequate to model the NRH 
instability. We derive the dispersion relation for the NRH 
instability in a simple case with the following assumptions. 
The CR are mono-energetic with momentum p. The zeroth 
order CR current jo, the wavenumber k and the zeroth order 
magnetic field Bo are all parallel. The background plasma is 
at rest to zeroth order, uo = 0. The first order perturbations 
to the magnetic field Bi, CR current ji, and plasma veloc- 
ity ui are all perpendicular to the zeroth order magnetic 
field and the wavevector k. Since the modes are transverse 
the plasma density is unperturbed to first order: pi = 0. 
The coupled linearised forms of equations (7), (13) and the 
Maxwell equation for dJi /dt are then: 

P^T = "Jo X Bi - ji x Bo + — (Bo.V)Bi 



dt 



/to 



gBi 

dt 



(Bo.V)ui 



dt 



c ec 6c 

--V x Gi - —B x ji Bi x j 

5 p p 



gGi 

dt 



cV x ji - vbGi 



(El) 



where Gi = (Aiv/3)ec J p 2 gidp. For circular polarisation 
any first order perturbation £i (Bi, ji, Gi or ui) satisfies 



dt 



-ojn x £i (n.V)£i = kn x £i (E2) 
where n is a unit vector in the direction of Bo, giving 



-j o B 1 +B j 1 + — Bi 
Mo 



-wji 



ck ecB . ecjo 
— n x Gi ji H Bi 

5 p p 



-wGi = cfcn x ji + i/jii x Gi 



(E3) 



As discussed in section 4 we set the scattering frequency vb 
equal to the CR Larmor frequency. For SNR conditions the 
NRH growth rate (~ T) is much less than the CR Larmor 
frequency, in which case the dispersion relation simplifies to 



, 2 2 
k V A 



a±T" 1 ± 



5i 
k 2 r g 



(E4) 



where va = / pfio is the Alfven speed, r g is the CR 

Larmor radius in the magnetic field Bo and T — WkBojo/p 
is the NRH growth rate where it dominates in the range 
r g l <C k <C F/va- <r± = ±1 according to the sense of the 
circular polarisation as determined by the sign of k. 

For wavelengths much shorter than the Larmor radius 
(kr g > 1) 



2 7 2 2 ti2 
UJ = rZ V A — CT±1 



(Eh) 



which represents a purely growing instability for wavenum- 
bers less than T/va in the appropriate polarisation a± = 1. 
Tension in the magnetic field correctly damps waves with 
wavenumbers greater than F/va (Bell 2004). The truncated 
tensor analysis is correct on scales smaller than the Larmor 
radius because CR trajectories are then relatively unaffected 
by perturbations in the magnetic field. This is the important 
regime in which the rapidly growing NRH instability ampli- 
fies the magnetic field. 

At wavelengths longer than the Larmor radius (kr g <C 
1) the approximate tensor expansion gives 



,2 2 , ■ ' 

= K Va ± i- 



(EG) 



In comparison the correct dispersion relation in this limit 
is ui 2 = k 2 v A — k 2 r 2 T 2 /5 for the resonant instability and 
uj 2 = k 2 VA + k 2 r 2 F 2 /5 for the non-resonant instability which 
in fact is stable for mono-energetic CR in this limit. The 
Alfven term k 2 v\ is negligible at long wavelengths so the 
tensor expansion gives a growing mode with 



i±l 



(kr B )r 



(E7) 



— uBi = kBoUi 



in both resonant and non-resonant polarisations. 

Figure 5 compares the approximate tensor dispersion 
relation (top row) with the correct dispersion relations of 
Bell (2004) (bottom row) for both the resonant instability 
(a± = — 1) that dominates for kr g < 1 and the rapidly grow- 
ing non-resonant NRH instability (a± = — 1) that dominates 
for kr g > 1. Crucially the tensor expansion accurately cal- 
culates the NRH growth rate in the range r" 1 < k < Y /va- 
Outside this range the growth rate is too small for the insta- 
bility to be effective. The tensor expansion reproduces the 
growth of the resonant instability for kr g < 1 although it 
incorrectly produces weak growth of the non-resonant insta- 
bility in this regime. The crucial point for the simulation 
is that the tensor expansion gives an accurate account of 
instability where the growth rate is large. 

Appendix F: Computational constraints 

The simulation makes heavy demands on computational 
resources because it models spatial scales encompassing the 
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wavelength of the fastest growing mode, the CR Larmor ra- 
dius, and the free propagation of CR ahead of the shock. It is 
three-dimensional in configuration space and models the CR 
distribution in momentum. The spatial cell-size Aa; must be 
small enough to allow the NRH instability to grow from its 
initially small scale (Ax < n/k ma x) and the timestep At 
must be short enough to resolve CR crossing one compu- 
tational cell (At < Ax/c). The simulation must be run for 
a time r at least 10 times the growth time of the fastest 
growing mode r = 107m„ a; and the length of the compu- 
tational grid Li I in the direction parallel to the shock nor- 
mal must be large enough to allow the CR to escape up- 
stream: Lii = 10cj max . Provided the boundary conditions 
are periodic in the directions perpendicular to the shock 
normal they can be much smaller than Lii but they must 
be able to accommodate a CR Larmor radius: L± = r g . 
The number of computational operations is proportional to 
N cornp = TL\\L 2 L /(AtAx 3 ), which from definitions and equa- 
tions presented above, is of order 



As discussed above the NRH instability grows strongly 
in a wavenumber range between r" 1 and k max . Saturation 
due to magnetic tension occurs on scalelengths comparable 
to the CR Larmor radius r g>sa t when B sa t/r g , S at = fJ-ojcn 
as discussed in section 10. The ratio of the saturated field 
B S at to the initial field B is B 3at /B = (2r g k max ) 1/2 where 
kmax and r g are defined in the initial field Bo- The equation 
for k max can be found in section 2, giving 



where the acceleration efficiency rj is also defined in section 
2. The Alfven Mach number Ma must be large to allow 
significant amplification of the magnetic field (B 3at 3> -Bo)- 
From equation (CI) large Ma imposes a heavy demand on 
computational resources, but from equation (C2) the cost 
can be minimised by making the shock velocity u s a large 
fraction of the speed of light. We initiate the simulation with 
an unrealistically large magnetic field of to reduce the 

Alfven Mach number from its typical value of ~ 1000 for 
shocks in young SNR. 




(F2) 
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